---
title: "7-3 all cause mortality"
description: |
Analyses of all cause mortality.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
toc-depth: 5
editor_options:
chunk_output_type: console
---
```{r}
#| label: pkgs
#| code-summary: Load packages
library (ASCOTr)
library (tidyverse)
library (lubridate)
library (kableExtra)
library (patchwork)
library (cmdstanr)
library (posterior)
library (bayesplot)
library (ggdist)
library (lme4)
library (broom)
library (broom.mixed)
library (bayestestR)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
all_dat <- readRDS (file.path (ASCOT_DATA, "all_data.rds" ))
# FAS-ITT
fas_itt_dat <- ASCOTr::: make_fas_itt_set (all_dat)
fas_itt_nona_dat <- fas_itt_dat |>
filter (! is.na (D28_death))
# ACS-ITT
acs_itt_dat <- ASCOTr::: make_acs_itt_set (all_dat)
acs_itt_nona_dat <- acs_itt_dat |>
filter (! is.na (D28_death))
# AVS-ITT
avs_itt_dat <- ASCOTr::: make_avs_itt_set (all_dat)
avs_itt_nona_dat <- avs_itt_dat |>
filter (! is.na (D28_death))
```
```{r}
#| label: stan-models
#| code-summary: Load models
logistic_mod <- compile_cmdstanr_mod (
file.path ("binary" , "logistic" ), dir = "stan" )
logisitc_site <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site" ), dir = "stan" )
logistic_site_epoch <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site_epoch" ), dir = "stan" )
```
```{r}
make_D28_death_table_A <- function (dat, format = "html" ) {
tab <- dat %>%
mutate (AAssignment = factor (
AAssignment,
levels = paste0 ("A" , 0 : 2 ),
labels = intervention_labels ()$ AAssignment)
) %>%
group_by (AAssignment) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
` Died within 28 days ` = sprintf ("%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
) %>%
bind_rows (
dat %>%
group_by (AAssignment = "Overall" ) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
` Died within 28 days ` = sprintf ("%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
)
) %>%
mutate (AAssignment = fct_inorder (AAssignment)) %>%
gather (key, value, - AAssignment, factor_key = T) %>%
spread (AAssignment, value)
colnames (tab)[1 ] <- "n ( \\ %)"
if (format == "latex" ) {
colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )
}
kable (tab,
format = format,
align = "lrrrrr" ,
escape = FALSE ,
booktabs = TRUE
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position" )
}
make_D28_death_table_C <- function (dat, format = "html" ) {
tab <- dat %>%
mutate (CAssignment = factor (
CAssignment,
levels = paste0 ("C" , 0 : 4 ),
labels = intervention_labels ()$ CAssignment)
) %>%
group_by (CAssignment) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
` Died within 28 days ` = sprintf ("%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
` Died within 28 days ` = sprintf ("%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
)
) %>%
mutate (CAssignment = fct_inorder (CAssignment)) %>%
gather (key, value, - CAssignment, factor_key = T) %>%
spread (CAssignment, value)
colnames (tab)[1 ] <- "n ( \\ %)"
if (format == "latex" ) {
colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )
}
kable (tab,
format = format,
align = "lrrrrrr" ,
escape = FALSE ,
booktabs = TRUE
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position" )
}
```
# Descriptive
## Anticoagulation
```{r}
#| label: tbl-domainC-po
#| code-summary: Relationship anti-coagulation to outcome
#| tbl-cap: Primary outcome by anti-coagulation intervention.
make_D28_death_table_C (
acs_itt_dat
)
save_tex_table (make_D28_death_table_C (
acs_itt_dat,
"latex" ),
"outcomes/secondary/7-3-anticoagulation-summary" )
```
## Antiviral
```{r}
#| label: tbl-domainA-po
#| code-summary: Relationship anti-viral to outcome
#| tbl-cap: Primary outcome by anti-viral intervention.
make_D28_death_table_A (
avs_itt_dat
)
save_tex_table (make_D28_death_table_A (
avs_itt_dat,
"latex" ),
"outcomes/secondary/7-3-antiviral-summary" )
```
## Age
```{r}
#| label: fig-age-po
#| code-summary: Relationship age to outcome
#| fig-cap: |
#| Relationship (logistic regression linear in age)
#| between age at entry and the primary outcome.
agedat <- fas_itt_dat %>%
dplyr:: count (D28_death, AgeAtEntry) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 0 ` + ` 1 ` + ` <NA> ` ,
p = ` 1 ` / (` 1 ` + ` 0 ` ))
agemod <- glm (
cbind (` 1 ` , ` 0 ` ) ~ AgeAtEntry,
data = agedat,
family = binomial ())
agedat <- agedat %>%
mutate (
ypred = predict (agemod, newdata = agedat, type = "response" )
)
p1 <- ggplot (agedat, aes (AgeAtEntry, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" )
p2 <- ggplot (agedat, aes (AgeAtEntry, p)) +
geom_point () +
geom_vline (xintercept = 60 , linetype = 2 ) +
geom_line (aes (y = ypred)) +
labs (y = "Proportion \n died by day 28" , x = "Age at entry" )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-age.pdf" ), p, height = 2 , width = 6 )
p
```
## Sex
```{r}
#| label: fig-sex-7-3
#| code-summary: Relationship sex to outcome
#| fig-cap: Mortality by sex.
tdat <- all_dat %>%
filter_fas_itt () %>%
dplyr:: count (Sex, D28_death) %>%
group_by (Sex) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
)
p1 <- ggplot (tdat, aes (Sex, n)) +
geom_col () +
labs (
y = "Number of \n participants" )
p2 <- ggplot (tdat, aes (Sex, p_1)) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ) +
ylim (0 , NA )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-sex.pdf" ), p, height = 2 , width = 6 )
p
```
## Oxygen
```{r}
#| label: fig-oxygen-7-3
#| code-summary: Relationship oxygen to outcome
#| fig-cap: Mortality by oxygen.
tdat <- all_dat %>%
filter_fas_itt () %>%
dplyr:: count (
supp_oxy2 = factor (supp_oxy2, labels = c ("Not required" , "Required" )),
D28_death) %>%
group_by (supp_oxy2) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
)
p1 <- ggplot (tdat, aes (supp_oxy2, n)) +
geom_col () +
labs (
y = "Number of \n participants" , x = "Supplemental oxygen" )
p2 <- ggplot (tdat, aes (supp_oxy2, p_1)) +
geom_point () +
labs (
y = "Proportion \n died by day 28" , x = "Supplemental oxygen" ) +
ylim (0 , NA )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-oxygen.pdf" ), p, height = 2 , width = 6 )
p
```
## Country
```{r}
#| label: fig-country-7-3
#| code-summary: Relationship country to outcome
#| fig-cap: Mortality by country.
tdat <- all_dat %>%
filter_fas_itt () %>%
dplyr:: count (Country = fct_infreq (PT_CountryName), D28_death) %>%
group_by (Country) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
)
p1 <- ggplot (tdat, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (tdat, aes (Country, p_1)) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Country of enrolment" ) +
ylim (0 , NA )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-country.pdf" ), p, height = 2 , width = 6 )
p
```
## Site
```{r}
#| label: fig-site-7-3
#| code-summary: Relationship site to outcome
#| fig-cap: Day 28 mortality by site within country.
tdat <- all_dat %>%
filter_fas_itt () %>%
dplyr:: count (
Country_lab = Country,
Site_lab = fct_infreq (Location),
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" )),
Site = PT_LocationName,
D28_death) %>%
group_by (Country, Site) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
) %>%
ungroup ()
p1 <- ggplot (tdat, aes (Site_lab, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (tdat, aes (Site_lab, p_1)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Site of enrolment" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p <- p1 / p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-cal-po
#| code-summary: Relationship calendar date to outcome
#| fig-cap: |
#| Relationship between calendar date and the primary outcome.
caldat <- all_dat %>%
filter_fas_itt () %>%
dplyr:: count (D28_death, yr = year (RandDate), mth = month (RandDate)) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (p = ` 1 ` / (` 1 ` + ` 0 ` ),
n = ` 1 ` + ` 0 ` + ` <NA> ` )
p1 <- ggplot (caldat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 ) +
ylim (0 , 0.15 )
p2 <- ggplot (caldat, aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p <- p2 | p1
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
# Analyses
## FAS-ITT
```{r}
res <- fit_primary_model (fas_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death" )
res$ drws$ OR <- res$ drws$ OR[- 1 ]
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-fas-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), FAS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-3-primary-model-fas-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-epoch-site-terms-fas-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-fas-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
::: {.callout-caution collapse="true"}
### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
::: {.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
### Posterior predictive
```{r}
#| label: primary-acs-itt-ppc
#| fig-cap: Posterior predictive distribution versus observed proportion (red diamond).
y_ppc <- res$ drws$ y_ppc
ppc_dat <- bind_cols (fas_itt_nona_dat, tibble (y_ppc = y_ppc))
grp_ppc <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
mean_y = mean (D28_death),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
}
plot_grp_ppc <- function (dat, lab = "" ) {
dat %>%
ggplot (aes (y = grp, xdist = rvar_mean_y_ppc)) +
stat_interval (size = 2 ) +
geom_point (aes (x = mean_y, y = 1 : nrow (dat)), colour = "red" , shape = 23 ) +
labs (
x = "Posterior predictive \n proportion" ,
y = lab)
}
ppc_A <- grp_ppc (sym ("AAssignment" ))
ppc_C <- grp_ppc (sym ("CAssignment" ))
ppc_ctry <- grp_ppc (sym ("Country" ))
ppc_epoch <- grp_ppc (sym ("epoch" ))
ppc_site <- ppc_dat %>%
group_by (Country = Country, grp = site) %>%
summarise (
mean_y = mean (D28_death),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
p0 <- plot_grp_ppc (ppc_A, "Antiviral \n intervention" ) + labs (x = "" )
p1 <- plot_grp_ppc (ppc_C, "Anticoagulation \n intervention" ) + labs (x = "" )
p2 <- plot_grp_ppc (ppc_ctry, "Country" ) + labs (x = "" )
p3 <- plot_grp_ppc (ppc_epoch, "Epoch" ) + labs (x = "" )
p4 <- plot_grp_ppc (ppc_site %>% filter (Country == "IN" ), "Sites \n India" ) + labs (x = "" )
p5 <- plot_grp_ppc (ppc_site %>% filter (Country == "AU" ), "Sites \n Australia" ) + labs (x = "" )
p6 <- plot_grp_ppc (ppc_site %>% filter (Country == "NP" ), "Sites \n Nepal" )
p7 <- plot_grp_ppc (ppc_site %>% filter (Country == "NZ" ), "Sites \n NZ" )
p <- ((p3 | p0 / p1 / p2) /
( (p4 | p5) / (p6 | p7) +
plot_layout (heights = c (3 , 1 )) ) ) +
plot_layout (heights = c (1 , 1.5 ), guides = "collect" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,
"7-3-primary-model-fas-itt-ppc.pdf" )
ggsave (pth, p, width = 6 , height = 5.5 )
p
```
## ACS-ITT
```{r}
#| label: fit-model-acs-itt
#| code-summary: Fit primary model
res <- fit_primary_model (acs_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death" )
res$ drws$ OR <- res$ drws$ OR[- 1 ]
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-acs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), ACS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-3-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-epoch-site-terms-acs-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
::: {.callout-caution collapse="true"}
### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
::: {.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
## AVS-ITT
### Pre-specified
- excludes epoch
```{r}
res <- fit_primary_model (
avs_itt_nona_dat,
logisitc_site,
ctr = contr.equalprior,
outcome = "D28_death" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" , "ctry" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 1 )
)
res$ drws$ OR <- res$ drws$ OR[- 1 ]
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Oxygen requirement" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-pre-spec
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-3-primary-model-avs-itt-summary-table-pre-spec" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
::: {.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
::: {.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["tau_site" ])
```
:::
### Reduced Model
- excludes epoch, region, and site
```{r}
res <- fit_primary_model (
avs_itt_nona_dat,
logistic_mod,
ctr = contr.equalprior,
outcome = "D28_death" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 )
)
res$ drws$ OR <- res$ drws$ OR[- 1 ]
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Oxygen requirement" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" )
```
```{r}
#| eval: false
XX <- res$ dat$ X
XX_m <- colMeans (XX)
XX <- cbind (XX[, 1 : 6 ], sweep (XX[, 7 : 12 ], 2 , XX_m[7 : 12 ]))
mdat <- res$ dat
mdat$ X <- XX
snk <- capture.output (
mfit <- logistic_mod[["sample" ]](
data = mdat,
seed = 1234 ,
adapt_delta = 0.95 ,
refresh = 0 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
chains = 4 ,
parallel_chains = min (4 , parallel:: detectCores ())
)
)
mpars <- mfit$ metadata ()$ model_params
keep <- mpars[! grepl ("(_raw|epsilon_)" , mpars)]
mdrws <- as_draws_rvars (mfit$ draws (keep))
names (mdrws$ beta) <- colnames (mdat$ X)
mdrws$ OR <- exp (mdrws$ beta[! grepl ("(Intercept|rand)" , names (mdrws$ beta))])
median (mdrws$ beta)
```
```{r}
#| label: odds-ratio-summary-table-avs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-3-primary-model-avs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-avs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
::: {.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
::: {.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
# End of script.